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We introduce a time-dependent projected Gross-Pitaevskii equation to describe a partially con- 
densed homogeneous Bose gas, and find that this equation will evolve randomised initial wave 
functions to equilibrium. We compare our numerical data to the predictions of a gapless, second 
order theory of Bose-Einstein condensation [S. A. Morgan, J. Phys. B 33, 3847 (2000)], and find 
that we can determine a temperature when the theory is valid. As the Gross-Pitaevskii equation 
is non-perturbative, we expect that it can describe the correct thermal behaviour of a Bose gas as 
long as all relevant modes are highly occupied. Our method could be applied to other boson fields. 



The achievement of Bose-Einstein condensation (BEC) 
in a dilute gas offers the possibility of studying the dy- 
namics of a quantum field at finite temperatures in the 
laboratory [0J2). However, direct numerical simulation 
of the full equations of motion for such systems is well 
beyond the capability of today's computers. Even equi- 
librium calculations in the region of a phase transition 
require non-perturbative methods, meaning that fully 
quantal treatments are unfeasible. 

At finite temperature, when there are an appreciable 
number of non-condensed particles, the fully quantal sec- 
ond order theory of Morgan || (and other equivalent 
treatments should be sufficient for an accurate de- 

scription of many properties of the dilute Bose gas in 
equilibrium and away from the region of critical fluctu- 
ations. Dynamical treatments are much harder, and in 
general require significant approximations. For example, 
calculations have been performed for small systems ||; 
with a restricted number of modes [Q; and for the dy- 
namics of condensate formation where the ground state 
is assumed to grow adiabatically || . 

The Gross-Pitaevskii equation (GPE) has been used 
to predict the properties of condensates near T = 0, 
when there are very few non-condensate atoms present. 
Both statically and dynamically it has shown excellent 
agreement with experiment ]9|-[lH. It has been argued, 
however, that the GPE can be used to describe the dy- 
namics of a Bose-Einstein condensate at finite tempera- 
ture flp^-prj. In the limit where the modes of the sys- 
tem are highly occupied (Nk ^> 1), the classical fluctu- 
ations of the field overwhelm the quantum fluctuations, 
and these modes may therefore be represented by a co- 
herent wave function. This is analogous to the situation 
in laser physics, where the highly occupied laser modes 
can be well described by classical equations. 

Using this argument, Damle et al. have performed 
calculations of the approach to equilibrium of a near 
ideal superfluid [ fl5| , and similar approximations to other 
quantum field equations have been successful elsewhere 
p6[ . References |l^-|l9|] also use the GPE to represent the 



classical modes of a Bose-condensed system. The main 
advantage of this method is that realistic calculations, 
while still a major computational issue, are feasible — 
methods for solving the GPE are well developed. Also, 
as the GPE is non-perturbative it should be possible to 
study the region of the phase transition, as long as the 
condition Nk ^ 1 is satisfied. 

There are, however, problems associated with the 
GPE. It is a classical equation, and so in equilibrium it 
will satisfy the equipartition theorem — all modes of the 
system will contain an energy fc^T. Thus, if we couple a 
system to a heat bath, and solve the equation with infi- 
nite accuracy, we will observe an ultra-violet catastrophe. 
Also, the higher the energy of any given mode, the lower 
its occupation will be in equilibrium, and eventually the 
criterion Nk 3> 1 will no longer be satisfied. For these 
low occupation modes a form of kinetic equation is more 
appropriate. The solution to both of these problems is to 
introduce a cutoff in the modes represented by the GPE. 

Our theoretical approach begins with the full operator 
equation for the Bose field with two-body interactions 
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where Uq = Airna/m is the effective interaction strength 
at low momenta, a is the s-wave scattering length, and 
to is the particle mass. The route to the usual GPE is 
to assume that the full field operator can be replaced by 
a wave function tp(r) — i.e. that all quantum fluctuations 
can be neglected. We proceed instead by defining a pro- 
jection operator V such that 
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where the region C is determined by the requirement 
that (a^ak) lj and the set {<^k} defines some basis 
in which the field operator is approximately diagonal at 
the boundary of C. For these modes, the quantum fluctu- 
ation part of the projected field operator can be ignored, 
and so we replace <2k — > Ck and write 



1 



1>(r) = Ck0k(r). (3) 
kec 

Defining the operator Q = 1 — V and Q$>(v) = fj(r), 
operating on Eq. (jl|) with V and taking the mean value 
results in what we call the finite temperature GPE 

^^P=ff ^(r) + ^{|^(r)| 2 V(r)} 

+ U Q f {2|VXr)| 2 (r)(r)) + ^(r) 2 (r}t( r ))} 

+ U T {V*(r)(7)(r)j7(r)) + 2^(r)<^(r)q(r)>} 

+ (7 P{(7}t(r)^(r)^(r)}}. (4) 

This describes the full dynamics of the region C and its 
coupling to an effective heat bath rj(r), which in prin- 
ciple can be described using a form of quantum kinetic 
theory. The finite temperature GPE is discussed in detail 
in Ref. §(J. 

In this letter, however, we wish to show that the GPE 
alone can describe the evolution of general configurations 
of the coherent region C towards an equilibrium that can 
be parameterised by a temperature. We therefore ignore 
all terms involving fj(r) in Eq. (gj) and concentrate on 
the first line, which we call the projected GPE. Although 
this equation is both unitary and reversible, we expect 
it to evolve general states to equilibrium, because de- 
terministic nonlinear systems exhibit chaotic, and hence 
ergodic, behaviour if more than a few degrees of free- 
dom are present plfl . This is confirmed by our numerical 
simulations and forms the main result of this letter. 

The projected GPE describes a microcanonical system. 
However, if the region C is large, then fluctuations in en- 
ergy and particle number in the grand canonical ensemble 
would be small. Hence we expect the final equilibrium 
state of the projected GPE to be similar to that of the 
finite temperature GPE coupled to a bath ?y(r) with the 
appropriate chemical potential and temperature. This 
does not affect the main result of this letter, however, 
which is simply that equilibrium is attained. The detailed 
non-equilibrium dynamics of the system will depend on 
the exchange of energy and particles between C and the 
bath, and this will be addressed in future work. 

We have performed simulations for a fully three- 
dimensional homogeneous Bose gas with periodic bound- 
ary conditions. This choice has been made to simplify 
the projection operation that must be carried out. In 
this case the single particle states are plane waves, and 
the effect of a condensate is simply to mix modes of mo- 
menta p and —p. This allows us to apply the projector 
cleanly in momentum space, which is easily accessible by 
fast fourier transform. In principle there is no barrier to 
performing the same computation in a trap — in practice 
the projection operation is much more time consuming. 

The dimensionless equation we compute is 

% dj^ L = _^ (f) + Cnl V\^(r)\ 2 ip(r), (5) 



where we have defined J <i 3 r|-0(r)| 2 = 1. The nonlin- 
ear constant is C n i = 2mNUo/h 2 L, where N is the to- 
tal number of particles in the volume, and L is the pe- 
riod of the system. Our dimensionless parameters are 
f = r/L, wave vector k = kL, energy e — e/el, and time 
t = ELt/h, with £l = h 2 /(2mL 2 ). 

The calculations presented here have been performed 
with C n i — 2000, and the projector V chosen such that 
all modes have |k| < 15 x 2n/L. This means that most 
of the states contained in the calculation are phonon- 
like for large condensate fraction. We note that while 
the number of states in the problem is fixed, the nonlin- 
ear constant only determines the ratio of NUq/L. This 
means that for a given value of C n i, we are free to choose 
N, Uq and L such that our condition |ck| 2 3> 1 is always 
satisfied for a given physical situation. In particular we 
can choose 87 Rb atoms with N — 5 x 10 5 and L » 17 /zm 
to give a number density of about 10 14 cm -3 — similar 
parameters to current experiments in traps. 
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FIG. 1. Condensate fraction plotted against total energy 
after each individual simulation has reached equilibrium. The 
barely discernable vertical lines on each point indicate the 
magnitude of the fluctuations. 

We begin our simulations with wave functions far from 
equilibrium with a chosen total energy E. They have a 
flat distribution in fc-space out to some maximum mo- 
mentum determined by E, and the phase of each mo- 
mentum component is chosen at random. These initial 
states are then evolved for a time period of r = 0.4, by 
which stage equilibrium appears to have been reached. 
We determine the properties of the system at equilibrium 
by assuming that the ergodic theorem applies, and time- 
averaging over 50 wave functions from the last St — 0.1 of 
the simulation. We find that the equilibrium properties 
depend only on the total energy — they are independent 
of the details of the initial wave function. 

Strong evidence that the simulations have reached 
equilibrium is given by the time dependence of the con- 
densate population. For all energies this settles down to 
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an average value that fluctuates by a small amount, and 
the results are presented in Fig. [j]. Further support is 
provided by the distribution of the particles in momen- 
tum space. Rather than using the plane-wave basis, we 
transform the wave functions into the quasiparticle ba- 
sis of quadratic Bogoliubov theory, the sole parameters 
of the transformation being the condensate fraction and 
the nonlinear constant C n \. We then average the popu- 
lations of the quasiparticles states over time and angle 
to produce a one dimensional plot, and the results are 
shown in Fig. |^. 
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FIG. 2. Plots of the equilibrium Bogoliubov quasiparticle 
distributions averaged over time and angle for four different 
total energies. Squares, E = 1600; crosses, E = 2000; circles, 
E — 3200; dots, E = 4600. The mean condensate occupation 
for all four distributions is off-axis. 

The GPE is the high occupation limit of the full equa- 
tion for the Bose field operator. Therefore, in equilib- 
rium we expect the mean occupation of the quasiparti- 
cle mode k to be the classical limit of the Bose-Einstein 
distribution i.e. the equipartition relation 

k B T 



(Nk) 



(6) 



Since we can determine the Bogoliubov occupation (Nk) 
from our simulation data, we can attempt to fit this dis- 
tribution to a dispersion relation for Sk , and hence deter- 
mine the temperature. 

In the limit of large condensate fraction (N )/N ~ 1, 
we expect the Bogoliubov dispersion relation to be a good 
estimate of the energies. The Bogoliubov transformation 
approximates the many-body Hamiltonian by a quadratic 
form, which can be diagonalised exactly. The eigenstates 
are quasiparticles, and in our dimensionless units the dis- 
persion relation takes the form 



e k = U 4 + 2C n 



N 



1/2 



(7) 



Manipulating Eq. (^J) and measuring the excitation spec- 
trum relative to the condensate we find 



f 



( N 



N 



IW (N ) 



(8) 



where f = k B T / (Nei) is our dimensionless temperature, 
and the second term on the RHS arises from the differ- 
ence between the condensate energy and the chemical 
potential of the system. By comparing the curve of this 
relation with that of Eq. (Q) , a temperature can be deter- 
mined. For the E = 1400 simulation we find f = 0.0284 
gives an excellent fit, and this is shown in Fig. ||. At 
higher simulation energies, however, the shape of Eq. (|^) 
no longer agrees with Eq. (Q) and we must use a more 
sophisticated theory to predict the dispersion relation. 
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FIG. 3. Comparison of dispersion relations for E — 1400 
with (No)/N = 0.912. The dots are a plot of the RHS of 
Eq. (Jsh , and the lines plot ik/T&f The solid line is for the 
Bogoliubov dispersion relation with Tat ~ 0.0284, while the 
dashed line is for the ideal gas with Tjj t = 0.0223. 

As the occupation of the quasiparticle modes becomes 
significant (in this case more than a few percent), the 
cubic and quartic terms of the many-body Hamiltonian 
that were neglected in the Bogoliubov transformation be- 
come important. In Ref. || Morgan develops a consis- 
tent extension of the Bogoliubov theory to higher order 
that leads to a gapless excitation spectrum. This theory 
treats the cubic and quartic terms of the Hamiltonian 
using perturbation theory in the quasiparticle basis. Ex- 
pressions for the energy-shifts of the excitations are given 
in Sec. 6.2 of Ref. ||, and we have calculated these shifts 
for our simulations, with typical results plotted in Fig. |^. 

The energy spectrum predicted by the second order 
theory for the E = 4000 simulation is in good agree- 
ment with the quasiparticle populations extracted from 
the simulations, and is a significantly better fit than the 
Bogoliubov theory of Eq. (0) . The validity of the second 
order theory is constrained by the requirement pi 



k B T 
n U 



(n a 3 ) 1/2 « 1, 



(9) 



3 



where no is the condensate density. For the results of 
Fig. [| with E — 4000, this parameter is 0.14 and so we 
are beginning to probe the boundary of validity of the 
theory. At higher E, the shifts it predicts at low k are 
of the order of the unperturbed energies, and the results 
become unreliable. In this region higher order terms are 
important and the second order theory can no longer be 
expected to give good results. 
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FIG. 4. Comparison of dispersion relations for E — 4000 
with (No)/N = 0.279. The dots are a plot of the RHS of 
Eq. (H) and the lines plot e^/faf The dashed line is for the 
Bogoliubov dispersion relation with Tfl t = 0.193, and the solid 
line is for the second order theory of Ref. J3j with fat = 0.201. 

In summary for the system with Cnl = 2000, Bogoli- 
ubov theory gives a good prediction of the energy spec- 
trum for simulations with total energies E < 1600, while 
the predictions of second order theory are good up until 
about E « 4000. We would like to point out, however, 
that as the GPE is non-perturbative we expect it will be 
valid up to and beyond the transition region as long as 
the condition Nk ^ 1 is satisfied. 

In addition to the results described above, we have also 
run simulations with C n \ = 10000 and carried out an 
identical analysis. We have found that the results from 
evolving the GPE are qualitatively the same, and for very 
large condensate fractions Bogoliubov theory accurately 
predicts the energy spectrum accurately. However, it ap- 
pears that the second order theory develops a gap in the 
energy spectrum in systems with a momentum cutoff. 
This feature is yet to be understood. 

In this Letter we have presented results for some of the 
equilibrium properties of the homogeneous gas. Other 
properties such as fluctuations, coherence lengths, as well 
as the non-equilibrium dynamics will be considered else- 
where. We would like to emphasise that this method re- 
lies on the lowest energy modes of the system being clas- 
sical in nature, and thus cannot handle situations where 
strong quantum fluctuations are important. 



In conclusion, we have presented evidence that the pro- 
jected Gross-Pitaevskii equation is a good approximation 
to the classical modes of a Bose gas. We have described 
how to carry out the projection technique in the homoge- 
neous case with periodic boundary conditions, and have 
shown that starting with a randomised wave function 
with a given energy, the projected GPE evolves towards 
an equilibrium state. We have analysed the numerical 
data in terms of the gapless, finite temperature theory 
of Ref. H in the classical limit, and found that both the 
occupation and energies of the quasiparticles agree quan- 
titatively with the predictions. 

We would like to thank R. J. Ballagh and C. W. Gar- 
diner for useful discussions. This work was financially 
supported by St John's College and Trinity College, Ox- 
ford, and the UK-EPSRC. 



[1 
[2 
[3 
[4 

[5 

[6. 

ir 

[8 
[9 
[10 

[11 

[12 
[13 

[i4; 
[is; 

[16 
[17 
[18 

[19 
[20 



M. Anderson et al. Science 269, 198 (1995). 

K. B. Davis et al. Phys. Rev. Lett. 75, 3969 (1995). 

S. A. Morgan, J. Phys. B 33, 3847 (2000). 

P. O. Fedichev and G. V. Shlyapnikov, Phys. Rev. A 58, 

3146 (1998). 

S. Giorgini, Phys. Rev. A 61, 063615 (2000). 

P. D. Drummond and J. F. Corney, Phys. Rev. A. 60, 

R2661 (1999). 

M. Holland, K. Burnett, C. W. Gardiner, J. I. Cirac and 

P. Zoller, Phys. Rev. A 54, R1757 (1996). 

C. W. Gardiner, M. D. Lee, R. J. Ballagh, M. J. Davis 

and P. Zoller, Phys. Rev. Lett. 81, 5266 (1999). 

M. Edwards and K. Burnett, Phys. Rev. A 51, 1382 

(1995). 

P. A. Ruprecht, M. J. Holland, K. Burnett and M. Ed- 
wards, Phys. Rev. A 51, 4704 (1995). 

F. Dalfovo, S. Giorgini, L. P. Pitaevskii and S. Stringari, 
Rev. Mod. Phys. 71, 463 (1999). 

B. V. Svistunov, J. Moscow Phys. Soc. 1, 373 (1991). 
Yu. Kagan, B. V. Svistunov, and G. V. Shlyapnikov, Zh. 
Eksp. Teor. Fiz. 101, 528 (1992) [Sov. Phys. JETP 75, 
387 (1992)]; Yu. Kagan and B. V. Svistunov, Phys. Rev. 
Lett. 79, 3331 (1997). 

R. J. Marshall, G. H. C. New, K. Burnett, and S. Choi, 
Phys. Rev. A 59, 2085 (1999). 

K. Damle, S. N. Majumdar, and S. Sachdev, Phys. Rev. 
A 54, 5037 (1996). 

G. D. Moore and N. Turok, Phys. Rev. D 55, 6538 (199 7). 
M. J. Bijlsma and H. T. C. Stoof, |cond-mat/0007026| . 
K. Goral, M. Gajda, and K. Rzazewski, Opt. Express 8, 
92 (2001). 

A. Sinatra, C. Lobo, and Y. Castin, |cond-mat/0101210[ 
M. J. Davis, R. J. Ballagh, and K. Burnett, cond 



mat/0107515, to appear in J. Phys. B. 
[21] L. E. Reichl, A Modern Course in Statistical Physics, 
(University of Texas Press, Austin, 1980). 



4 



